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Investigation of Allan Variance for Determining Noise Spectral Forms 
With Application to Microwave Radiometry 

by 

William D. Stanley 
Abstract 

An investigation of the Allan Variance method as a possible means for 
characterizing fluctuations in radiometric noise diodes has been performed. The goal is to 
separate fluctuation components into white noise, flicker noise, and random-walk noise. 
The primary means is by discrete-time processing, and the study focused primarily on the 
digital processes involved. Noise satisfying the requirements was generated by direct 
convolution. Fast Fourier transformation (FFT) processing in the time domain, and FFT 
processing in the frequency domain. Some of the numerous results obtained are presented 
along with the programs used in the study. 

Introduction 

The application of the Allan Variance method to determine the stability of precision 
oscillator circuits has been well documented in the literature (refs. 1—7). R. W. Lawrence 
of NASA Langley Research Center (ref. 8) has proposed to apply the same concept to the 
determination of the long-term stability of noise diodes for precision radiometric 
applications. This concept has been investigated in this report and some of the findings 
have been documented here. 



Sample Mean and Variance of Random Process 

Assume a random process x(t ) defined over some time interval, and let x k (t) 
represent the time average of this function over the interval t k < t < t k + r as illustrated in 

figure 1(a). This operation can be performed by the integrate-and-dump process, which 
can be described as 

x k (t) = ^j' t+T x(t)dt (1) 


The value x k (x) is the time-average or mean over an interval with duration r. 

Assume now that the function is sampled at equally-spaced intervals T seconds 
apart as shown in figure 1(b). For a given interval having N discrete samples, a discrete- 
time approximation to the continuous-time formulation of (1) is the sample mean as defined 
by 

i k+N - 1 

( 2 ) 

This operation can also be referred to in signal-processing terminology as a sum-and-dump 
algorithm. 

The sample variance of this process is denoted as v k (N,x) and is given by 


. i *+a/-i 

v * (Ar,T)= ivM 


(3a) 


N - 1 


*+ AM 


fk+N - 1 


S I 


\ * 


(3b) 


The proceeding formulation is based on the statistics of a so-called sample 
population, from which information about the process can be determined from the sample 
mean and sample variance. The term "Allan Variance" refers to a specific manner in which 
the sample variance is determined as a function of the integration time f to sort out 
processes exhibiting different types of power spectra. 
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Assume that the power spectrum S(f) is of the form 


S{f) = for 0</<oo (4) 

where K x is a constant. The theoretical work of Allan and others establishes the basis for 
the variance v(N, t) to have the form 

v(N,T) = K 2 r k ~ l (5) 


where K 2 is another constant. 

Some of the underlying theory, however, involves improper integrals and diverging 
functions so it is not always clear how the proceeding process will behave with finite 
summations and discrete-time processes. 


From (4) and (5), several forms were of primary interest in this investigation. They 


are as follows: 



White Noise: 

ii 

Q 

(6) 


v(N, t) = — 
T 

(7) 

Flicker Noise: 

S(f) = 

(8) 


v(N, X) = K 2 

(9) 

Random Walk Noise: 

S<f) = jL 

(10) 
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By determining v(N, t) as a function of t, it may be possible to sort out the spectral terms 
according to the behavior of the variance. 

It can be shown that the expected value of v k (N, x) is given by 

E[v k (N, t)] = x 2 - E[x k x t ] = x 2 — Vi x ( r) (12) 

where x 2 is the mean-square value of the process and 9l,(r) = is the 

autocorrelation function corresponding to a time shift of T. 

In general, the mean-square value can be expressed as 

+ (13) 

where cP is the variance and fi x is the mean or dc value. Similarly, the autocorrelation 
function can be expressed as 


9t,(T) = J?,(T) + AI* (14) 

where R x { r) is the autocovariance function. Substitution of (13) and (14) in (12) results 
in 


E[v k (N, r)] = a 2 , + - R x ( r) - v 2 

= o 2 x -R x (t) 


(15a) 

(15b) 


Thus, the effect of any dc component in the random process cancels out and does not affect 
the variance. 


For any realistic processes, R x (t) approaches zero as t increases. This means 
that the preceding expectation tends to approach 

E[v*(iV,T)]«^ (16) 

as t increases. 

After studying the pattern of the algorithm for different values of N , it was decided 
to restrict most of the study to a two-point variance as suggested by Lawrence. While it is 
possible to produce slightly lower fluctuations in the variance estimates by using more 
points in the algorithm, there is a tradeoff with respect to resolution. Further, the 
complexity of implementation of the algorithm increases markedly with increasing N. 
Thus, N = 2 was judged to be a good choice, and in subsequent expressions, the variance 
rotation will be simplified to v(N, t) = v(2, t) = v(r), with N = 2 understood. (Later in 
the report, N will be used as the number of points in a record length.) 

White Noise Reduction by Integration 

It is instructive to initially review the procedure by which the variance of a white 
noise process can be reduced by integration or averaging. Consider white noise having a 
flat power spectrum as shown in figure 2. It has been assumed here that the noise has been 
sampled at intervals of T seconds apart and, hence, the sampled spectrum is defined over 
the frequency range -1/2 T <f<\/2T. The power spectrum has a constant value 
S(f) = N 0 /2 over the region, and the sampling process creates aliases at integer multiples 
of the sampling frequency. The variance of the noise can be expressed as 

c 2 =— . 2B= N 0 B = — (17) 

x 2 2 T 
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Assume that the function x(t) is successively integrated or averaged over an interval x 
with x increasing. The equivalent noise bandwidth of the process is then 0.5/ x. Letting 
v(t) represent the variance of the filtered process, it can be readily shown that 

v(T) = -<7? (18) 

X 

The conclusion is that the variance of a white-noise process is inversely proportional to the 
integration time for an integrate and dump or sum and dump filter. The process is 
illustrated in figure 3. 


Two-Point Allan Variance 

Consider a random process x(t) for which N equally-spaced samples are taken 
at intervals of T seconds over a total epoch length of t p seconds; i.e., 


t p =NT 


(19) 


Assume that N is selected as an integer power of 2; i.e.. 


N= 2‘ 


( 20 ) 


where / is an integer. Consider the first two points x x and x 2 and let 


_ jc, + x, 

X = _1 L 


( 21 ) 


The sample variance v, , of this two-point distribution is 


-*) 2 = -x? 

* i i 


( 22 ) 


n :■ r 


Substitution of (21) in (22) and simplification lead to 



( 23 ) 


The preceding process is repeated at N! 2-1 intervals of width 2 T to yield a 
sequence of two-point variance samples as follows: 


1 , 

V 1.2 = ^- C ^ 3 ~^ 4 > 

(24) 

V l.3=^(*5-* 6 ) 2 

(25) 

V lj=^( X 2j-\ ~ X 2j) 2 

(26) 

1 ' .2 

v i.w/2 = 2^v-i -a: v) 

(27) 


Next, the mean value v, of the Nil variances is formed as follows: 

j Nl 2 


(28) 


The quantity v, is the average or mean variance associated with the short integration 
interval r = T. This constitutes the completion of the first cycle. 

The second cycle of the computational process is initiated by first doubling the 
"integration" or averaging interval to 2 T by successive linear averages. Let 


x,+x 2 
■* 2.1 - ~ 


(29) 
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(30) 


■*2,2 “ 


*3+ *4 


„ _ *2,-1+ *2,,- 

^ - 2 


_ *^-1 +X N 

x 2,sn ~ 2 


The variance estimates for the longer intervals are now computed by 


V 2 ,l =-(*2,1 *2,2 ) 2 


V 2 .2=-j(*2,3-*2,4) 2 


V 2 ,;=^(* 2 , 2 ;- l -* 2 , 2;) 2 


V 2,N/4 — 2 (*2,iV/2-l *2,M/2 ) 


The mean variance v 2 is now computed as 


N / 4 

V,=— £v 2 , y 


wt 


(31) 

(32) 

(33) 

(34) 

(35) 

(36) 

(37) 


Without showing the intermediate steps, there will be a total of L = log 2 N cycles. 
On the last cycle, the interval is NT/2 and there are only two samples in which the variance 
can be computed. Hence 


V L~ 2^ Xl,X 


(38) 
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Algorithm for Two-Point Variance 

The analysis in this section will be limited to the two-point variance algorithm, for 
which most of the data obtained are based. The basic process in cyclical form is illustrated 
in figure 4. Basically, there are L = log 2 N cycles with three nearly identical steps except 
for the first and last cycles. 

The first cycle consists of the following three steps: 

(1) transfer of N samples of x(t) to y(t ) preserve original array intact 

(2) computation of N/2 separate estimates of the variance v 1>; - by the two-point 

algorithm and 

(3) computation of the mean variance v from the N/2 separate variances based 
on an integration time T. 

The first step in the second though the L — 1 cycles consists of replacing the N 
points of y with N/2 points based on doubling the sampling time interval. Steps (2) and 

(3) are identical to those of the first cycle except that the number of points used in the 
separate variance and mean variance estimates become successively smaller in each cycle. 

The first two steps in the Lth or last cycle follow the preceding pattern except that 
there is only one variance estimate v L in this case. Thus, there is no mean variance 
compilation in this case and v L -v L . The equivalent integration time is NT / 2. 

As a further aid in visualizing the process, refer to the chart of figure 5 based on 
N = 16 points. Obviously, this record length would be far too short to obtain any 
statistical reliability, but the chart provides an interesting visual representation of the 
process. In this case L = log 2 16 = 4. The various integration times in this process are 
T, 2 T, 4T, and 8 T. There are four cycles and a total of 12 steps. 

A computer flow chart for the Allan Variance algorithm is shown in figure 6. Some 

of the notation varies slightly from that previously utilized due to the programming 
constraints. In fact, the two-dimensional variables y t J and v i } can be replaced by a one 

dimensional variable by replacing the array with new data following each cycle. 
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Statistical Analysis of the Allan Variance 


The extent to which the Allan variance algorithm is statistically significant will now 
be investigated. Let v, y represent the ith cycle of the algorithm based on the y'th segment, 

which will be denoted for this development simply as v(N, t) . The expression for v(N, t) 
is 

v(V,t) = ~— X(;y* - 50 2 (39) 

N - 1 i 

Assume that y(t) is gaussian distributed white noise with mean fi y and variance a y . 
The Allan variance can then be expressed as 

2 

v(w,t ) =-^a;L, (40) 


where Xn-i is a chi-squared distribution with N - 1 degrees of freedom as defined by 


*V-l=*l 2 +*2 + 


+ x 


N 


(41) 


Each x n is gaussian distributed with mean 0 and variance 1. 

The mean value and variance of the chi-square distribution are denoted, 
respectively, as \i , and o 2 2 , and they can be shown to be 

X tf — 1 X N — 1 



The expected value of v(N, r) is then given by 


E[v(N, r)] = /i v = — z- x (IV - 1) = c) 
N - 1 7 


(44) 


The variance of v(N, t) is 


E{[v(N, t)-fi v ] 2 ) 


2(<7 y 2 ) 2 (jV-l) 2g* 

(N- 1) 2 N - 1 


(45) 


The mean value v(t) based on N i variance computations is 


v(r) = ~^v k (N,T) 

i 


(46) 


It can be shown that 


E[V(T)]=a 2 y =n v (47) 

and 

2<j 4 

E[(v-h v ) 2 ] = (J^ = 2 — (48) 

v v Ni ( N _y) 

The preceding results have been analyzed to determine the various statistical 
parameters inherent in the Allan variance algorithm. A tabular formulation of these results 
is presented in table 1 for a record length of N points satisfying the requirement that 
N = Np, where N p is the number of points in each cycle chosen as an integer power of 2 
and L is the number of cycles in the algorithm. 

From the results of the table, it is clear that the mean value of the variance as well as 
the mean value of the averaged variance both decrease inversely with the order of the cycle. 
This represents the expected trend for white noise since the order of the cycle is 
proportional to the integration time. 
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The variance a * of the variance computation decreases inversely with the number 
of points in each cycle and with the square of the cycle number. The variance cr? of the 
averaged variance decreases inversely with the number of points in each cycle and with the 
order of the cycle. 

A statistical merit factor that was found to be very useful is the ratio <r 7 / fJ.-, i.e., 
the ratio of the standard deviation of the averaged variance to the value of the averaged 
variance. This parameter provides a measure of the relative uncertainty of the averaged 
variance. This parameter varies inversely with the square root of the number of points in a 
given cycle. In all cases, this parameter is smaller for small n, but it increases as the 
number of points used in the estimate decreases. In fact for a two-point algorithm, the 
merit factor reaches a level of V2 for the last cycle. This extreme case means that the 
relative uncertainty on the last cycle is about ±140 % of the estimated variance! 

A significant conclusion of the preceding analysis is that the record length must be 
considerably longer than the desired integration time over which the estimate is to be made. 
Depending on the desired accuracy, it could be required to be 100 or more times the desired 
integration time. Another way of looking at this is to say that the relative accuracy of the 
last portion of the record length is considerably uncertain. 

A similar tabulation for two-point variance is shown in table 2. To illustrate some 
actual numbers, table 3 was constructed for N = 64 and two-point variance. Finally, 
table 4 illustrates the results that would apply for four-point variance and for N = 64. 

Variance Analysis with Discrete Operations 

The original work performed by Allan and others primarily utilized Fourier 
theorems based on classical continuous functions. While this approach is more 
mathematically rigorous, it is intuitively void because of the divergence of some of the 
functions involved at f = 0 and/or / = 0. Examples of some of the pertinent Fourier 
pairs are shown in figure 7. Inasmuch as any modem measurements would likely be made 
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utilizing discrete sampling processes (possibly with an FFT), it seems prudent to 
investigate the phenomena involved based on a finite number of sampled points. 

Consider the three forms of noise spectra depicted in figure 8. The functions of (a), 
(b), and (c) are assumed, respectively, as power spectra of white noise, flicker noise, and 
random-walk noise. Each spectrum is presumed to be derived from a process sampled at a 
rate of f s = 1 / T. However, the power spectra fold about the folding frequency 
fo = fs /2 , and a unique spectrum exists only from dc to The effect of the Allan 

variance will be analyzed for each case based on a discrete summation of the power under 
each curve. 

White Noise. For white noise, the power spectrum will be expressed as 

S(f) = ~ ~ for 0 </</„ (49) 

J o 

where cr 2 is the variance of the process. Let v(r) represent the Allan Variance based on 
an integration time t, The corresponding frequency resolution A/", is 

V,=- (50) 

T, 


The variance of the unmodified function is determined by taking the area under the curve 
between 0 and f 0 and it is simply 


v(r 1 ) = ^-A/,+^-4f 1 + 

Jo Jo 


N h 2 A c Nh 2 

= — cr A/, = — —o 2 

f f 

do Jo 


'l' 
v T i y 


( 51 ) 


where N H - N / 2. In the initial case, N H A/, = f o and the result is a 2 as expected. 
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Next, assume that the integration time is doubled, i.e., let r, =2r,, and let 
Afj = 1 / t 2 = 1 / 2 t, . The corresponding frequency increment is now half as large as 
before, and the variance can now be written as 


V ( T 2) = ~J~ 4fl + + 

^ 0 Jo 


ml 

So 


A h = 


f ML' 


So 


2t, 


(52) 


It should be noted that f a is based on the original sampling rate and is a constant. Thus, 
the variance has been reduced by one-half and is cx 2 / 2. 

It is easy to generalize this process at this point and write for any r 

v(r) = — (53) 

T 


Thus, it is expected that the Allan variance is inversely proportional to the integration time 
for white noise. 

Ricker Noise . Now assume that the power spectrum is of the form 

SCO -4 o </</, (54) 

This form is depicted in figure (b). Because of the divergence of (45) at / = 0, the lowest 
frequency indicated is Af in which S(Af) = K/Af. 
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The variance in this case can be written as 


v( t) = S(Af)Af + 5(2 40 Af + L + S{NAf)Af 




1 1 


1 


A f 2Af 3A f 


+ .... 


N„Af 


w 


AT, | 1 + - + - + 
2 3 


+ ■ 


N, 


Hj 


= K X K 2 


(55) 


where K 2 is a constant representing the summation. The expected result is a constant 
value independent of the integration time! 

Random Walk Noise . Finally, assume that the spectrum is of the form 

S(/) = yf- 0 </</„ (56) 

This form is illustrated in figure (c). Again, because of the divergence at / = 0, the lowest 
frequency component is Af in which S(f) = K/(Af) 2 . 

The variance in this case can be written as 


v (t) = S(Af)Af + S(2Af)Af + S(3Af)Af + .... +S(NAf)Af 


= K, 


1 1 


[(A/ 2 ) (2A/f (3Af) : 


-+ .... + 


(W w A/J 


A f 


= K, 


,11 111 

4 9 N 2 \ Af 




A/ 


( 57 ) 
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This result indicates that the variance of random-walk noise should increase linearly with 
the integration time T. 

Resolution . It should be understood that as the integration time increases, the 
number of points used in the variance estimate decreases, and the resulting uncertainty of 
the estimate increases. A different way of looking at this is that the lowest frequency in the 
spectrum is A/“ = 1 / t and this is the region where the uncertainty is greatest. 
Consequently, the integration time must be much larger than the reciprocal of the lowest 
frequency at which an estimate is desired. 

Computer Programs Developed 

Numerous computer programs were written during the course of this study, but four 
particular programs deserve special attention due to their extensive utilization within the 
study. These four programs were written in QBASIC ("Quick Basic") and are identified by 
the following four titles: (1) ALLAN.BAS, (2) ALLCON.BAS, (3) ALLFFT.BAS, and 
(4) ALLFF1.BAS. The first one explicitly uses the name Allan in reference to Allan 
Variance, and the others use the abbreviation "ALL" for Allan at the beginning of the title. 
All of the programs utilize the Allan variance algorithm, but the last three provide much 
more. Each will be discussed briefly, but detailed program listings will be delayed until the 
last three sections. 

(1) ALLAN.BAS . This program provides a direct computation of the Allan Variance of a 
signal having 2 1 points, where / is an integer. In theory / could be set to any arbitrary 
integer, but a practical limit is imposed by the storage capacity of the computer used. The 
output is a tabulation of the Allan variance as a function of the integration time r. 
Although any input data could be utilized, most of the studies utilizing this program directly 
employed white noise. 

(2) ALLCON.BAS . This program performs a direct convolution in the time domain 
between white noise and a specified impulse response whose Fourier transform has the 
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proper spectral shape required to produce different types of noise. The key feature here is 
that all of the analysis is done directly in the time domain. The last part of the program 
tabulates the Allan variance as a function of the integration time. 

(3) ALLFFT.BAS . This program starts with an input random process and an impulse 
response as in the previous program but does not convolve in the time domain. Instead, 
the FFT of the input signal and the FFT of the impulse response are first computed. Next 
the FFT's are multiplied at each frequency using complex arithmetic, and the inverse FFT 
is computed. This signal is then applied to the Allan variance process, and the Allan 
variance as a function of the integration time is then computed. 

(4) ALLFFT1.BAS . This program starts with an input random signal and computes the 
FFT. Unlike the previous two programs, an impulse response is not provided. Instead, the 
FFT of the input signal is adjusted in the frequency domain to the desired spectral shape. 
The inverse FFT is then computed, and the Allan variance algorithm is applied. The output 
is a tabulation of the Allan variance as a function of the integration time. 

White Noise Simulation 

White noise can be directly simulated with either of the four programs, although it 
is on "overkill" with the latter three. The basis for the random signal in each case was the 
random number generator RND in QBASIC. This generator produces a random variable 
that is uniformly distributed from 0 to 1 with a mean value of 0.5 and a standard deviation 
of 0. 5 / V3 (or a variance of 1/12). After a reasonable number of averaging operations, a 
gaussian distribution may be assumed. In all simulations, the variable was shifted to 
produce a mean of zero and multiplied by 100. The result is a distribution having a mean 
value of zero and a variance of 2500/3 = 833.3. 

The experimental system employed in the laboratory utilized a sample time 
T = 0.2 s. Consequently, the output data format for integration time was adjusted to create 
values of time that were multiples of the sampling time. 
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A listing of the program ALLAN and the corresponding simulation results for white 
noise with N = 2 13 = 8192 points are shown in figure 9. An associated plot is shown in 
figure 10. For the initial integration time of 0.25, the computed Allan variance is about 
819.8, which differs from the theoretical value of 833.3 by about 1.62%. From the plot, 
the variance as a function of integration time follows closely the expected 1/r dependence 
up to about t = 100 s, at which time the variation begins to deviate significantly. 

Although "clumsy" because of the additional features, it was decided as a test to 
compare all four programs with white noise as the desired input variable. In this test, the 
number of points was limited to 2 10 — 1024 points due to memory restrictions. The 
program ALLAN was run again, and a program listing along with the output data are 
shown in figure 11. The corresponding plot is shown in figure 12. Once again, the 
derivation from the expected behavior becomes apparent as r approaches the limit of the 
time range. 

Next, the program ALLCON was adjusted for white noise. This was achieved by 
setting the impulse response simply to H( 0) = 1 and H(N) = 0 for N * 0. The program 
listing and output data are shown in figure 13. These results are in perfect agreement with 
those of the preceding program and, therefore, were not plotted. 

The third run of the same input signal was made with ALLFFT using the same 
noise generator as in the preceding run. Once again, the impulse response was established 
as a single point with non-zero value, i.e. H( 0) = 1. The program listing is shown in 
figure 14 and the output data are shown in figure 15. Once again, the results are the same. 

The fourth run of the same data was made with ALLFF1. In this case, only the 
input random signal was specified in the time domain. The approach in this case was to 
specify the spectral shaping. This is achieved with the use of the functions HM{M) and 
HP(M). The quantity HM(M) represents the frequency domain filter magnitude 
response and HP(M) represents the phase response. These functions were established 
for this run as follows: 
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HM(M) = 1 for 0<M < 2043 


HP{ 0) = 0 for 0 < M < 2043 (58) 

The program listing is shown in figure 16, and the data are shown in figure 17. 
The results are again in agreement. 

Flicker Noise Simulation 

The generation of flicker noise and the measurement of the Allan Variance was 
achieved separately with the three programs ALLCON, ALLFFT, and ALLFFT1. Flicker 
noise possesses a power spectrum proportional to 1 // . The process used to create the 
desired noise is based on the following Fourier transform pair: 



Since convolution in the time domain is equivalent to multiplication in the frequency 

domain, one has the choice of either convolving the white noise at the input with an 
impulse response proportional to 1 / Vt (or 1 / -Jn ) or by multiplying the white noise 

spectrum by a transfer function proportional to 1 / ■yff (or 1 / Vm). 

When the linear spectrum is modified by a function proportinal to 1 / -yff , the 

power spectrum becomes proportional to 1// , which is the desired outcome. 
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First, ALLCON was used for an 8192-point signal. The impulse response was 


defined as 


( 0 ) = 0 

H(N)=l/4N for 1 < N <8192 (60) 

Note that it is necessary to define H( 0) separately due to the divergence of the spectral 
form at N = 0 . What this means is that the impulse response is invalid at the first point, 
but this limitation did not seem to create any problems with the results. 

The program listing and a table of the Allan variance versus integration time are 
shown in figure 18. A corresponding plot is shown in figure 19. The ideal expected 
values of the Allan variance should be essentially a constant, and this trend is definitely 
observable up to about r = 25 s. For the last decade or so, however, the process is erratic 
as expected. 

Next, a run was taken with ALLCON based on a 1024-point array corresponding to 
r= 102.4 s. The program listing and data shown in figure 20. The corresponding plot is 
shown in figure 21. The nearly constant behavior of the variance at small values of x is 
evident, but the results become unstable near the end of the time interval. 

The preceding signal was then applied to the program ALLFFT. To avoid overlap, 
the signal was first padded with zeros. Specifically the following input was created: 

X(N) = random noise for 0 < N < 1023 

= 0 for 1024 <N< 2047 (61) 

H( 0) = 0 (62) 


20 


i 1 1 ' 



H{N) = \l4N for 1 < iV < 1023 


= 0 for 1024 <N < 2047 (63) 

A program listing is shown in figure 22, and the output data are shown in figure 
23. These results are in exact agreement with the direct convolution approach in the time 
domain. Note, however, that due to the padding, a 2048-point FFT was used to generate a 
1024-point result for variance analysis. 

Finally, ALLFFT1 was used for a 1024-point signal. The input random signal was 
padded with zeros from N = 1024 to 2047 as in the previous case. However, the spectral 
shaping in this case was performed in the frequency domain by the use of the following 
transfer function: 

HM( 0) = 0 (64) 

HM(M) = 1 / VM for 1 < M < 1023 (65) 

HM (1024) = 0 (66) 

HM(M) = HM( 2048 - M ) for 1025 < M < 2047 (67) 

and 

HP(M) = 0 for 0<M£ 2047 (68) 

Note that the strategy in the frequency domain is to create a function whose 
magnitude response is an even function about the folding frequency, which corresponds to 
M = 1024 . This is necessary for the corresponding impulse response to be a real function. 
If the phase response had been specified as anything other than zero, it would be necessary 
to force it to be an odd function about M = 1024 . 

The program listing is shown in figure 24, and the output data are shown in 

figure 25. In this case, the level of the output signal is quite different than with the use of 

ALLCON and ALLFFT. The reason is that there was no attempt made to adjust the 
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spectral level to correspond with that obtained in the preceding two runs. However, the 
relative level is the quantity of interest, and a plot is shown in figure 26. The relative level 
is observed to be nearly a constant over a wide range of % as expected. 

Random Walk Noise Simulation 

Continuing in the tradition established earlier, random walk noise was generated in 
the time domain by convolving white noise with a constant value for the response. A 
program utilizing ALLCON for a duration of 409.6 s and the output data are shown in 
figure 27. The plot of these data is shown in figure 28. The upward trend of the variance 
as r increases is evident. 

The array was reduced to 1024 points, constituting a time of 102.4 s, and the 
program and data are shown in figure 29. The corresponding plot is shown in figure 30. 

Next, the program ALLFFT was applied to the same random signal. In this case, a 
2048-point FFT was used with zero padding as discussed earlier. The output data are 
shown in figure 31. For all practical purposes, these results are in perfect agreement with 
those of the preceding case. 

Finally, ALLFFT 1 was used for the same signal by performing spectral alteration in 
the frequency domain. Random walk noise possesses a power spectrum proportional to 
l// 2 . The linear spectrum can then be adjusted by multiplying by 1// or, equivalently in 
form 1/A/. Specifically, the following transfer function was used: 


HM(0) = 0 

(69) 

HM( M) = l/M for 1<A/<1023 

(70) 

HM( 1024) = 0 

(71) 

HM(M) = HM (2048 - M) for 1025 < M <2047 

(72) 

22 



hi 8 1 



HP(M) = 0 for 0 < M < 2047 


(73) 


A program listing is shown in figure 32, and the output data are shown in figure 
33. A plot of the results is shown in figure 34. As in the case of flicker noise, the level is 
quite different here but the relative shape displays the correct pattern for most of the time 
range. 
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Summary and Conclusions 


The Allan variance process has been investigated as a means for characterizing the 
power spectra of random noise processes. Extensive simulations were performed using 
white noise, flicker noise, and random walk noise. Discrete time operations utilizing both 
convolution and FFT processes were utilized, and several computer programs were written 
to facilitate the study. 

In general, the Allan variance estimates followed the theoretical processes expected 
over much of the integration time range. However, two conditions were clearly evident 
from the results: (1) the record length must be considerably longer than the interval over 
which the integration time results must be evaluated. Said differently, the last 90 percent or 
so of the record length exhibits wide fluctuations and (2) the record length must be 
considerable larger than the reciprocal of the lowest frequency at which the spectrum is to 
be evaluated. 
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Table 1. Statistical properties as a function of general parameters 
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n = variable assuming integer values N '~ 1 for 1 <1<L. 

N = total number of points chosen as an integer power of N p . 

N p = number of points in each cycle chosen as an integer power of 2. 
N d = N/N p = largest number of points used in computing v . 




Table 2. Statistical properties for general two-point algorithms 
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n = variable assuming integer values 2' -1 for \ <1<M. 

N = total number of points chosen as an integer power of 2. 

N p =2 = number of points in each variance computation. 

N d = largest number of points used in computing mean variance. 
= N/2 
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Table 3. Statistical properties for N = 64 and N p = 2 
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N = 64 , N p = 2 N d = 6 4/2 = 32. 


28 


71 









































Table 4. Statistical properties for N = 64 and N p = 4 
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Figure 2. White noise spectral form after sampling. 
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Integration Time 
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V|_,j = ^(Yi - Y2) 2 for j = 1 , 1 (or N/N) 
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Figure 4. Algorithm for two-point variance computations. 
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Figure 6. Flow chart for two-point variance computation. 
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CLS 

DIM X(8192 ) , Y (8192 ) , V(8192) 

PRINT "PROGRAM ALLAN TO DETERMINE SAMPLE VARIANCE": PRINT 
PRINT "L", "INT TIME", "VARIANCE": PRINT 
NP = 8192: LP = 13 
FOR N = 0 TO NP - 1 

X(N) = 100 * (RND - .5): Y(N) = X(N) 

NEXT N 

FOR N = NP TO 1 STEP -1 
Y(N) * Y(N - 1) 

NEXT N 
N1 « NP 

FOR L = 1 TO LP 
NX * N1 / 2 
FOR N = 1 TO N1 

V(N) = .5 * (Y( 2 * N - 1) - Y( 2 * N)) " 2 
NEXT N 
VBAR = V(l) 

FOR N = 2 TO N1 

VBAR = VBAR + (V(N) - VBAR) / N 
NEXT N 

FOR N s 1 TO N1 

Y(N) = .5 * ( Y( 2 * N - 1) + Y( 2 * N)) 

NEXT N 

PRINT L, .2 * (2 ~ (L - 1)), VBAR 

NEXT L 

BEEP 


PROGRAM ALLAN TO DETERMINE SAMPLE VARIANCE 


L INT TIME VARIANCE 


1 

.2 

819.8405 

2 

.4 

422.2106 

3 

.8 

214.7032 

4 

1.6 

103.3248 

5 

3.2 

52.65189 

6 

6.4 

27.02878 

7 

12.8 

14.13634 

8 

25.6 

6.824592 

9 

51.2 

2.61482 

10 

102.4 

1.387154 

11 

204.8 

.3692694 

12 

409.6 

.442937 

13 

819.2 

.603497 


Figure 9. Program and data for white noise with t = 819.2 s and ALLAN. 
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CLS 

DIM X(8192), Y (8192 ) , V(8192) 

PRINT "PROGRAM ALLAN TO DETERMINE SAMPLE VARIANCE": PRINT 
PRINT "L", "INT TIME", "VARIANCE": PRINT 
NP * 1024: LP * 10 
FOR N = 0 TO NP - 1 

X(N) = 100 * ( RND - .5): Y(N) = X(N) 

NEXT N 

FOR N = NP TO 1 STEP -1 
Y(N) = Y(N - 1) 

NEXT N 
N1 = NP 

FOR L = 1 TO LP 
N1 = N1 / 2 
FOR N = 1 TO N1 

V(N) = .5 * ( Y( 2 * N - 1) - Y(2 * N)) “ 2 
NEXT N 
VBAR = V(l) 

FOR N = 2 TO N1 

VBAR = VBAR + (V(N) - VBAR) / N 
NEXT N 

FOR N = 1 TO N1 

Y(N) = .5 * (Y(2 * N - 1) + Y(2 * N) ) 

NEXT N 

PRINT L, .2 * (2 “ (L - 1)), VBAR 

NEXT L 

BEEP 


PROGRAM ALLAN TO DETERMINE SAMPLE VARIANCE 
L INT TIME VARIANCE 


1 

.2 

761.5204 

2 

.4 

406.3706 

3 

.8 

212.9356 

4 

1.6 

98.75756 

5 

3.2 

44.2822 

6 

6.4 

13.76261 

7 

12.8 

18.59924 

8 

25.6 

8.317465 

9 

51.2 

3.744104 

10 

102.4 

.2659905 


Figure 1 1 . Program and data for white noise with x = 102.4 s and ALLAN. 
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CLS 

DIM X(8192), 11(8192), Y(8192), V(8192) 

PRINT "PROGRAM ALLCON USING DIRECT CONVOLUTION": PRINT 

PRINT "L", "INT TIME", "VARIANCE": PRINT 

NP = 1024: LP ■ 10 

FOR N « 0 TO NP - 1 

X(N) = 100 * (RND - .5) 

H(N) = 0: Y(N) = 0 
NEXT N 
H(0) = 1 

FOR N = 0 TO NP - 1 
FOR I = 0 TO N 

Y(N) = Y(N) + X(I) * H(N - I) 

NEXT I 
NEXT N 

FOR N = NP TO 1 STEP -1 
Y(N) = Y(N - 1) 

NEXT N: N1 = NP 

FOR L = 1 TO LP: N1 = N1 / 2 

FOR N = 1 TO N1 

V(N) = .5 * ( Y( 2 * N - 1) - Y( 2 * N) ) ~ 2 
NEXT N 
VBAR = V(l) 

FOR N = 2 TO N1 

VBAR = VBAR + (V(N) - VBAR) / N 
NEXT N 

FOR N ■ 1 TO N1 

Y(N) = .5 * ( Y( 2 * N - 1) + Y(2 * N) ) 

NEXT N 

PRINT L, .2 * (2 ~ (L - 1)), VBAR 
NEXT L: BEEP 


PROGRAM ALLCON USING DIRECT CONVOLUTION 


L INT TIME VARIANCE 


1 

.2 

761.5204 

2 

.4 

406.3706 

3 

.8 

212.9356 

4 

1.6 

98.75756 

5 

3.2 

44.2822 

6 

6.4 

13.76261 

7 

12.8 

18.59924 

8 

25.6 

8.317465 

9 

51.2 

3.744104 

10 

102.4 

.2659905 


Figure 13. Program and data for white noise with x = 102.4 s and ALLCON. 
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CLS 

DIM X(2048) , X1R(2048) , X1I(2048), X2R(2048), X2I(2048) 

DIM H(2048) , H1R(2048), H1I(2048), H2R(2048), H2I(2048) 

DIM Y(2048) , Y1R(2048), Y1I(2048), Y2R(2048), Y2I(2048), V(2048) 

PRINT "PROGRAM ALLFFT USING FFT TO PERFORM CONVOLUTION": PRINT 

PRINT "N", "INT TIME", "VARIANCE": PRINT : T = . 2 

NP = 2048: LP » 11: NH = NP / 2: PI « 3.14159: U = 2 * PI / NP 

H(0) = 1: FOR N = 1 TO NH - 1: H(N) = 0: NEXT N 

FOR N = 0 TO NH - 1: X(N) = 100 * (RND - .5): NEXT N 

FOR N * NH TO NP - 1: X(N) = 0: H(N) ■ 0: NEXT N 

FOR M « 0 TO NP - 1: X1R(M) = X(M): XII (M) = 0: NEXT M: IA = NP / 2: 

FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U): ZI = -SIN(IC * U) 

FOR M = IC TO ID - 1 

AR = X1R(M + IC): AI * X1I(M + IC) 

BR = ZR * X1R(M + ID) - ZI * X1I(M + ID) 

BI ■ ZI * X1R(M + ID) + ZR * XII (M + ID) 

X2R(M) * AR + BR: X2I(M) = AI + BI 

X2R(M + NH) » AR - BR: X2I(M + NH) = AI - BI: NEXT M 

IC = ID: ID * ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: X1R(M) = X2R(M): XII (M) = X2I(M): NEXT M: NEXT 

FOR M = 0 TO NP - 1: H1R(M) = H(M) : H1I(M) = 0: NEXT M: IA = NP / 2: 

FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = -SIN(IC * U) 

FOR M = IC TO ID - 1: AR = H1R(M + IC): AI = H1I(M + IC) 

BR = ZR * H1R(M + ID) - ZI * H1I(M + ID) 

BI = ZI * H1R(M + ID) + ZR * H1I(M + ID) 

H2R(M) = AR + BR: H2I(M) = AI + BI 

H2R(M + NH) = AR - BR: H2I(M + NH) * AI - BI : NEXT M 

IC = ID: ID ■ ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: H1R(M) = H2R(M) : H1I(M) = H2I(M): NEXT M: NEXT 

FOR M = 0 TO NP - 1 

Y1R(M) = X1R(M) * H1R(M) - X1I(M) * H1I(M) 

Y1I(M) = X1R(M) * H1I(M) + XII (M) * HlR(M): NEXT M: IA = NP / 2: IB 
FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = SIN(IC * U) 

FOR N = IC TO ID - 1: AR = Y1R(N + IC): AI = Y1I(N + IC) 

BR = ZR * Y1R(N + ID) - ZI * Y1I(N + ID) 

BI = ZI * Y1R(N + ID) + ZR * Y1I(N + ID) 

Y2R(N) = AR + BR: Y2I(N) = AI + BI 

Y2R(N + NH) = AR - BR: Y2I(N + NH) = AI - BI : NEXT N 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR N = 0 TO NP - 1: Y1R(N) = Y2R(N) : Y1I(N) = Y2I(N): NEXT N: NEXT 

FOR N = NP TO 1 STEP -1: Y(N) = Y1R(N - 1) / NP: NEXT N: N1 = NH 

FOR L = 1 TO LP - 1: Nl = Nl / 2 

FOR N = 1 TO Nl: V(N) = .5 * (Y(2 * N - 1) - Y(2 * N)) ~ 2: NEXT N: 

VBAR = V(l) 

FOR N = 2 TO Nl: VBAR = VBAR + (V(N) - VBAR) / N: NEXT N 

FOR N = 1 TO Nl: Y(N) = .5 * (Y(2 * N - 1) + Y(2 * N)): NEXT N 

PRINT L, T * (2 “ (L - 1)), VBAR: NEXT L: BEEP 


IB = 1 
L 

IB = 1 
L 

= 1 


Figure 14. Program and data for white noise with x = 102.4 s and ALLFFT. 
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PROGRAM ALLFFT USING FFT 

N INT TIME 

.2 
.4 
.8 
1.6 

3.2 
6.4 
12.8 
25.6 

51.2 

10 102.4 


TO PERFORM CONVOLUTION 

VARIANCE 

761.5205 

406.3707 

212.9357 

98.75756 

44.28219 

13.76262 

18.59924 

8.317466 

3.744104 

.2659903 


Figure 15. Data for white noise with x = 102.4 s and ALLFFT. 
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CLS 

DIM X(2048), X1R( 2048 ) , X1I(2048), X2R(2048) / X2I(2048), HM(2048), HP(2048) 
DIM Y(2048) , Y1R(2048) , Y1I(2048) / Y2R(2048), Y2I(2048), V(2048) 

PRINT "PROGRAM ALLFFT1 USING FFT WITH SPECTRUM SHAPING": PRINT 

PRINT "N", "INT TIME", "VARIANCE": PRINT : T = .2 

NP = 2048: LP = 11: NH = NP / 2: PI * 3.14159: U = 2 * PI / NP 

FOR N = 0 TO NH - 1: X(N) = 100 * (RND - .5): NEXT N 

FOR N a* NH TO NP - 1: X(N) = 0: NEXT N 

FOR M = 1 TO NH - 1: HM(M) = 1: HP(M) = 0: NEXT M 

HM( 0 ) = 1: HP ( 0 ) = 0: HM(NH) = 1: HP ( NH ) = 0 

FOR M = NH + 1 TO NP - 1: HM(M) = HM(NP - M) : HP(M) = 0: NEXT M 

FOR M = 0 TO NP - 1: X1R(M) = X(M) : X1I(M) = 0: NEXT M: IA = NP / 2 : IB = 1 

FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = -SIN(IC * U) 

FOR M = IC TO ID - 1 

AR = X1R(M + IC): AI = X1I(M + IC) 

BR = ZR * X1R(M + ID) - ZI * X1I(M + ID) 

BI - ZI * X1R(M + ID) + ZR * XII (M + ID) 

X2R(M) = AR + BR: X2I(M) = AI + BI 

X2R(M + NH) = AR - BR: X2I(M + NH) = AI - BI: NEXT M 

IC = ID: ID * ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: X1R(M) = X2R(M): XII (M) = X2I(M): NEXT M: NEXT L 

FOR M = 0 TO NP - 1: HR = HM(M) * COS(HP(M)): HI = HM(M) * SIN(HP(M)) 

Y1R(M) = X1R(M) * HR - X1I(M) * HI 

Y1I(M) = X1R(M) * HI + XII (M) * HR: NEXT M: IA = NP / 2 : IB = 1 
FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = SIN(IC * U) 

FOR N = IC TO ID - 1: AR = Y1R(N + IC): AI = Y1I(N + IC) 

BR = ZR * Y1R(N + ID) - ZI * Y1I(N + ID) 

BI = ZI * Y1R(N + ID) + ZR * Y1I(N + ID) 

Y2R(N) ■ AR + BR : Y2I(N) = AI + BI 

Y2R(N + NH) = AR - BR: Y2I(N + NH) = AI - BI : NEXT N 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR N = 0 TO NP - 1: Y1R(N) = Y2R(N) : Y1I(N) = Y2I(N): NEXT N: NEXT L 

FOR N = NP TO 1 STEP -1: Y(N) = Y1R(N - 1) / NP: NEXT N: N1 = NH 

FOR L = 1 TO LP - 1 : Nl = Nl / 2 

FOR N = 1 TO Nl: V(N) = .5 * (Y(2 * N - 1) - Y(2 * N)) ~ 2: NEXT N: 

VBAR = V(l) 

FOR N = 2 TO Nl: VBAR = VBAR + (V(N) - VBAR) / N: NEXT N 

FOR N = 1 TO Nl: Y(N) = .5 * (Y(2 * N - 1) + Y(2 * N)): NEXT N 

PRINT L, T * (2 ~ (L - 1)), VBAR: NEXT L: BEEP 


Figure 16. Program for white noise with x = 102.4 s and ALLFFT 1 . 
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PROGRAM ALLFFT1 USING FFT WITH SPECTRUM SHAPING 


N INT TIME VARIANCE 


10 


.2 

761.5205 

.4 

406.3707 

.8 

212.9357 

1.6 

98.75756 

3.2 

44.28219 

6.4 

13.76262 

12.8 

18.59924 

25.6 

8.317466 

51.2 

3.744104 

102.4 

.2659903 


Figure 17. Data for white noise with x = 102.4 s and ALLFFT1 . 
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CLS 

DIM X( 8192 ) , H( 8192 ) , Y(8192), V(8192) 

PRINT "PROGRAM ALLCON USING DIRECT CONVOLUTION": PRINT 
PRINT "L", "INT TIME", "VARIANCE": PRINT 
NP = 4096: LP = 12 
FOR N = 0 TO NP - 1 

X(N) = 100 * (RND - .5): Y(N) = 0: NEXT N 

FOR N = 1 TO NP - 1: H(N) = 1 / SQR(N): NEXT N: H(0) = 0 

FOR N = 0 TO NP - 1 

FOR I = 0 TO N 

Y(N) = Y(N) + X(I) * H(N - I) 

NEXT I 
NEXT N 

FOR N = NP TO 1 STEP -1 
Y(N) = Y(N - 1) 

NEXT N: N1 = NP 

FOR L = 1 TO LP: N1 = N1 / 2 

FOR N = 1 TO N1 

V(N) = .5 * (Y(2 * N - 1) - Y( 2 * N)) ** 2: NEXT N 
VBAR = V(l) 

FOR N = 2 TO N1 

VBAR = VBAR + (V(N) - VBAR) / N: NEXT N 
FOR N = 1 TO N1 

Y(N) = .5 * ( Y( 2 * N - 1) + Y( 2 * N)): NEXT N 
PRINT L, .2 * (2 * (L - 1)), VBAR 
NEXT L: BEEP 


PROGRAM 

ALLCON USING DIRECT 

CONVOLUTION 

L 

INT TIME 

VARIANCE 

1 

.2 

465.0879 

2 

.4 

483.4081 

3 

.8 

511.6693 

4 

1.6 

599.207 

5 

3.2 

747.6578 

6 

6.4 

800.0687 

7 

12.8 

1009.559 

8 

25.6 

1217.807 

9 

51.2 

179.5898 

10 

102.4 

87.83582 

11 

204.8 

956.2097 

12 

409.6 

3456.044 


Figure 18. Program and data for flicker noise with x = 409.2 s and ALLCON. 
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CLS 

DIM X(8192), 11(8192), Y(8192), V(8192) 

PRINT "PROGRAM ALLCON USING DIRECT CONVOLUTION": PRINT 
PRINT "L", "INT TIME", "VARIANCE": PRINT 
NP = 1024: LP = 10 
FOR N = 0 TO NP - 1 

X(N) = 100 * (RND - .5): Y(N) = 0: NEXT N 

FOR N = 1 TO NP - 1: H(N) = 1 / SQR(N) : NEXT N: H(0) = 0 

FOR N = 0 TO NP - 1 

FOR I = 0 TO N 

Y(N) = Y(N) + X( I ) * H(N - I) 

NEXT I 
NEXT N 

FOR N = NP TO 1 STEP -1 
Y(N) = Y(N - 1) 

NEXT N: N1 = NP 

FOR L * 1 TO LP: N1 = N1 / 2 

FOR N = 1 TO N1 

V(N) = .5 * ( Y( 2 * N - 1) - Y( 2 * N)) “ 2 
NEXT N 
VBAR = V(l) 

FOR N = 2 TO N1 

VBAR = VBAR + (V(N) - VBAR) / N 
NEXT N 

FOR N = 1 TO N1 

Y(N) = .5 * ( Y( 2 * N - 1) + Y( 2 * N)) 

NEXT N 

PRINT L, .2 * (2 * (L - 1)), VBAR 
NEXT L: BEEP 


PROGRAM ALLCON USING DIRECT CONVOLUTION 


INT TIME VARIANCE 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 


.2 

437.8973 

.4 

440.3997 

.8 

532.6196 

1.6 

517.8219 

3.2 

800.6263 

6.4 

846.5082 

12.8 

1505.267 

25.6 

1135.48 

51.2 

5.159271 

102.4 

60.71339 


Figure 20. Program and data for flicker noise with x = 102.4 s and ALLCON. 
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CLS 

DIM X(2048 ) , X1R(2048) , X1I(2048), X2R(2048), X2I(2048) 

DIM H( 2048 ) , H1R( 2048) , H1I(2048) / H2R(2048), H2I(2048) 

DIM Y( 2048 ) , Y1R(2048) , Y1I(2048), Y2R(2048), Y2I(2048), V(2048) 

PRINT "PROGRAM ALLFFT USING FFT TO PERFORM CONVOLUTION": PRINT 

PRINT "N", "INT TIME", "VARIANCE": PRINT : T = .2 

NP = 2048: LP = 11: NH = NP / 2: PI = 3.14159: U = 2 * PI / NP 

H( 0) = 0: FOR N = 1 TO NH - 1: H(N) = 1 / SQR(N) : NEXT N 

FOR N = 0 TO NH - 1: X(N) = 100 * (RND - .5): NEXT N 

FOR N = NH TO NP - 1: X(N) = 0: H(N) = 0: NEXT N 

FOR M = 0 TO NP - 1: X1R(M) = X(M): X1I(M) = 0: NEXT M: IA ■ NP / 2: 

FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = -SIN(IC * U) 

FOR M = IC TO ID - 1 

AR = X1R(M + IC): AI * XII (M + IC) 

BR ■ ZR * X1R(M + ID) - ZI * X1I(M + ID) 

BI = ZI * X1R(M + ID) + ZR * XII (M + ID) 

X2R(M) = AR + BR: X2I(M) = AI + BI 

X2R(M + NH) = AR - BR: X2I(M + NH) = AI - BI : NEXT M 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: X1R(M) = X2R(M): X1I(M) = X2I(M): NEXT M: NEXT 

FOR M = 0 TO NP - 1: H1R(M) = H(M) : H1I(M) = 0: NEXT M: IA = NP / 2: 

FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U): ZI = -SIN(IC * U) 

FOR M = IC TO ID - 1: AR = H1R(M + IC): AI = H1I(M + IC) 

BR = ZR * H1R(M + ID) - ZI * H1I(M + ID) 

BI = ZI * H1R(M + ID) + ZR * H1I(M + ID) 

H2R(M) = AR + BR: H2I(M) = AI + BI 

H2R(M + NH) = AR - BR: H2I(M + NH) = AI - BI: NEXT M 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: H1R(M) = H2R(M) : H1I(M) = H2I(M): NEXT M: NEXT 

FOR M = 0 TO NP - 1 

Y1R(M) = X1R(M) * H1R(M) - X1I(M) * H1I(M) 

Y1I(M) = X1R(M) * H1I(M) + XII (M) * H1R(M): NEXT M: IA = NP / 2: IB 
FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = SIN(IC * U) 

FOR N = IC TO ID - 1: AR = YiR(N + IC): AI = Y1I(N + IC) 

BR = ZR * Y1R(N + ID) - ZI * Y1I(N + ID) 

BI = ZI * Y1R(N + ID) + ZR * Y1I(N + ID) 

Y2R(N) = AR + BR: Y2I(N) = AI + BI 

Y2R(N + NH) = AR - BR: Y2I(N + NH) = AI - BI: NEXT N 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR N = 0 TO NP - 1: Y1R(N) = Y2R(N) : Y1I(N) = Y2I(N): NEXT N: NEXT 

FOR N = NP TO 1 STEP -1: Y(N) = Y1R(N - 1) / NP: NEXT N: N1 = NH 

FOR L « 1 TO LP - 1: N1 = N1 / 2 

FOR N = 1 TO Nl: V(N) = .5 * (Y(2 * N - 1) - Y(2 * N)) “ 2: NEXT N: 

VBAR = V(l) 

FOR N = 2 TO Nl: VBAR = VBAR + (V(N) - VBAR) / N: NEXT N 

FOR N = 1 TO Nl: Y(N) = .5 * (Y(2 * N - 1) + Y(2 * N)): NEXT N 

PRINT L, T * (2 “ (L - 1)), VBAR: NEXT L: BEEP 


IB = 1 
L 

IB = 1 
L 

= 1 


Figure 22. Program for flicker noise with x = 102.4 s and ALLFFT. 
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PROGRAM ALLFFT USING FFT TO PERFORM CONVOLUTION 
N INT TIME VARIANCE 


.2 

437.8973 

.4 

440.4008 

.8 

532.6213 

1.6 

517.8229 

3.2 

800.6278 

6.4 

846.5118 

12.8 

1505.262 

25.6 

1135.482 

51.2 

5.160076 

102.4 

60.71655 


Figure 23. Data for flicker noise with x = 1 02.4 s and ALLFFT. 
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CLS 

DIM X( 2048) , X1R( 2048 ) , X1I(2048), X2R(2048), X2l(2048), HM(2048), HP(2048) 
DIM Y ( 2048 ) , Y1R(2048) , Y1I(2048), Y2R(2048), Y2I(2048), V(2048) 

PRINT "PROGRAM ALLFFT1 USING FFT WITH SPECTRUM SHAPING": PRINT 

PRINT "N", "INT TIME", "VARIANCE": PRINT : T = .2 

NP = 2048: LP = 11 : NH = NP / 2 : PI = 3.14159: U = 2 * PI / NP 

FOR N = 0 TO NH - 1: X(N) = 100 * (RND - .5): NEXT N 

FOR N = NH TO NP - 1: X(N) = 0: NEXT N 

FOR M = 1 TO NH - 1: HM(M) = 1 / SQR(M) : HP(M) = 0: NEXT M 

HM( 0 ) = 0: HP(0 ) = 0: HM(NH) = 0: HP(NH) = 0 

FOR M = NH + 1 TO NP - 1: HM(M) = HM(NP - M): HP(M) = 0: NEXT M 

FOR M = 0 TO NP - 1: X1R(M) = X(M) : XlI(M) = 0: NEXT M: IA = NP / 2 : IB = 1 

FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = -SIN(IC * U) 

FOR M = IC TO ID - 1 

AR = X1R(M + IC): AI = X1I(M + IC) 

BR = ZR * X1R(M + ID) - ZI * X1I(M + ID) 

BI = ZI * X1R(M + ID) + ZR * XII (M + ID) 

X2R(M) = AR + BR: X2I(M) = AI + BI 

X2R(M + NH) = AR - BR: X2I(M + NH) = AI - BI : NEXT M 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: X1R(M) = X2R(M): X1I(M) = X2I(M): NEXT M: NEXT L 

FOR M = 0 TO NP - 1: HR = HM(M) * COS(HP(M)): HI = HM(M) * SIN(HP(M)) 

Y1R(M) = X1R(M) * HR - X1I(M) * HI 

Y1I(M) = X1R(M) * HI + X1I(M) * HR: NEXT M: IA = NP / 2: IB = 1 

FOR L = 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U): ZI = SIN(IC * U) 

FOR N = IC TO ID - 1: AR = Y1R(N + IC): AI = Y1I(N + IC) 

BR = ZR * Y1R(N + ID) - ZI * Y1I(N + ID) 

BI = ZI * Y1R(N + ID) + ZR * Y1I(N + ID) 

Y2R(N) = AR + BR: Y2I(N) = AI + BI 

Y2R(N + NH) = AR - BR: Y2I(N + NH) = AI - BI : NEXT N 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR N = 0 TO NP - 1: Y1R(N) = Y2R(N) : Y1I(N) = Y2I(N): NEXT N: NEXT L 

FOR N = NP TO 1 STEP -1: Y(N) = Y1R(N - 1) / NP: NEXT N: N1 = NH 

FOR L = 1 TO LP - 1: N1 = N1 / 2 

FOR N = 1 TO Nl: V(N) = .5 * (Y(2 * N - 1) - Y(2 * N) ) * 2: NEXT N: 

VBAR = V(l) 

FOR N = 2 TO Nl: VBAR = VBAR + (V(N) - VBAR) / N: NEXT N 

FOR N = 1 TO Nl: Y(N) = .5 * (Y(2 * N - 1) + Y(2 * N)): NEXT N 

PRINT L, T * (2 “ (L - 1) ), VBAR: NEXT L: BEEP 


Figure 24. Program for flicker noise with x = 102.4 s and ALLFFT1 . 
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PROGRAM ALLFFT1 USING FFT WITH SPECTRUM SHAPING 


N 


INT TIME VARIANCE 


10 


.2 

.4 

.8 

1.6 

3.2 

6.4 
12.8 
25.6 

51.2 

102.4 


1.219962 

1.214046 

1.171432 

1.17254 

.9797807 

.7463165 

1.43384 

1.68266 

.9090664 

.3821753 


Figure 25. Data for flicker noise with x = 102.4 s and ALLFFT1 . 
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CLS 

DIM X(8192), H(8192), ¥(8192), V(8192) 

PRINT "PROGRAM ALLCON USING DIRECT CONVOLUTION": PRINT 
PRINT "L", "INT TIME", "VARIANCE": PRINT 
NP = 4096: LP = 12 
FOR N = 0 TO NP - 1 

X(N) = 100 * (RND - .5): Y(N) = 0: H(N) = 1: NEXT N 
FOR N = 0 TO NP - 1 
FOR I - 0 TO N 

Y(N) = Y(N) + X( I ) * H(N - I) 

NEXT I 
NEXT N 

FOR N s NP TO 1 STEP -1 
Y(N) = Y(N - 1) 

NEXT N: N1 = NP 

FOR L = 1 TO LP: N1 = N1 / 2 

FOR N = 1 TO N1 

V(N) = .5 * ( Y( 2 * N - 1) - Y( 2 * N) ) “ 2 
NEXT N 
VBAR = V(l) 

FOR N = 2 TO N1 

VBAR = VBAR + (V(N) - VBAR) / N 
NEXT N 

FOR N = 1 TO N1 

Y(N) = .5 * ( Y( 2 * N - 1) + Y(2 * N)) 

NEXT N 

PRINT L, .2 * (2 ~ (L - 1)), VBAR 
NEXT L: BEEP 


PROGRAM ALLCON USING DIRECT CONVOLUTION 


INT TIME VARIANCE 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 


.2 

.4 

.8 

1.6 

3.2 

6.4 

12.8 

25.6 

51.2 

102.4 
204.8 

409.6 


404.4132 

621.1478 

1031.323 

1914.698 

3831.717 

8453.081 

18015.49 

25252.21 

28800.24 

90826.44 

309915.6 

699118.1 


Figure 27. Program and data for radom walk noise with x = 409.6 s and ALLCON. 
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CLS 

DIM X( 1024 ) , H( 1024 ) , Y(1024), V(1024) 

PRINT "M", "INT TIME", "VARIANCE" 

PRINT 

NP = 1024 

FOR N * 0 TO 511 

X(N) = 100 * ( RND - .5): H(N) = 1: Y(N) = 0 
NEXT N 

FOR N = 512 TO NP - 1 

X(N) = 100 * (RND - .5): H(N) = 1: Y(N) = 0 
NEXT N 

FOR N = 0 TO NP - 1 
FOR I = 0 TO N 

Y(N) = Y(N) + X(I) * H(N - I) 

NEXT I 
NEXT N 

FOR N = NP TO 1 STEP -1 
Y(N) = Y(N - 1) 

NEXT N 
N1 = NP 

FOR M = 1 TO 10 
N1 = N1 / 2 
FOR N = 1 TO N1 

V(N) = .5 * ( Y( 2 * N - 1) - Y( 2 * N) ) “ 2 
NEXT N 
VBAR = V(l) 

FOR N a 2 TO N1 

VBAR = VBAR + (V(N) - VBAR) / N 
NEXT N 

FOR N = 1 TO N1 

Y(N) ■ .5 * (Y( 2 * N - 1) + Y(2 * N) ) 

NEXT N 

PRINT M, .2* (2 * (M- 1)), VBAR 

NEXT M 

BEEP 


INT TIME VARIANCE 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 


.2 

394.7494 

.4 

610.4443 

.8 

1093.376 

1.6 

1712.925 

3.2 

3951.782 

6.4 

12355.76 

12.8 

15130.16 

25.6 

38779.27 

51.2 

64480.73 

102.4 

117919.4 


Figure 29. Program and data for random walk noise with x = 102.4 s and ALLCON. 
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Figure 30. Plot for results of Figure 29. 


CLS • i. 

DIM X(2048) , X1R(2048) , XII(2048), X2R(2048), X2I(2048) 

DIM H(2048), H1R( 2048) , H1I(2048), H2R(2048), H2I(2048) 

DIM Y(2048), Y1R(2048 ) , ¥11(2048), Y2R(2048), Y2I(2048), V(2048) 
PRINT "PROGRAM ALLFFT USING FFT TO PERFORM CONVOLUTION": PRINT 
PRINT "N", "INT TIME", "VARIANCE": PRINT : T = .2 
NP » 2048: LP - 11: NH - NP / 2: PI * 3.14159: U = 2 * PI / NP 
FOR N ■ 0 TO NH - 1: X(N) ■ 100 * (RND - .5): H(N) * 1: NEXT N 
FOR N s NH TO NP - 1: X(N) = 0: H(N) ■ 0: NEXT N 

FOR M = 0 TO NP - 1: X1R(M) * X(M) : XII (M) = 0: NEXT M: IA ■ NP / 2: 

FOR L * 1 TO LP: IC = 0: ID - IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI = -SIN(IC * U) 

FOR M * IC TO ID - 1 

AR * X1R(M + IC): AT = XII (M + IC) 

BR * ZR * X1R(M + ID) - ZI * X1I(M + ID) 

BI ■ ZI * X1R(M + ID) + ZR * XII (M + ID) 

X2R(M) - AR + BR: X2I(M) » AI + BI 

X2R(M + NH) = AR - BR: X2I(M + NH) » AI - BI : NEXT M 

IC » ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: X1R(M) = X2R(M): X1I(M) = X2I(M): NEXT M: NEXT 

FOR M = 0 TO NP - 1: H1R(M) = H(M) : H1I(M) = 0: NEXT M: IA = NP / 2 : 

FOR L » 1 TO LP: IC = 0: ID = IA 

FOR K = 1 TO IB: ZR * COS(IC * U): ZI = -SIN(IC * U) 

FOR M = IC TO ID - 1: AR = HlRfM + IC): AI = HlI(M + IC) 

BR » ZR * H1R(M + ID) - ZI * H1I(M + ID) 

BI s ZI * H1R(M + ID) + ZR * H1I(M + ID) dJ 

H2R(M) = AR + BR: H2I(M) a AI + BI 

H2R(M + NH) * AR - BR: H2I(M + NH) = AI - BI: NEXT M 

IC = ID: ID » id + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M = 0 TO NP - 1: H1R(M) = H2R(M): H1I(M) = H2I(M): NEXT M: NEXT 

FOR M * 0 TO NP - 1 

Y1R(M) a X1R(M) * H1R(M) - XlI(M) * H1I(M) 

Y1I(M) a X1R(M) * H1I(M) + XII (M) * H1R(M) : NEXT M: IA = NP / 2 : IB 

FOR L » 1 TO LP: IC = 0: ID = IA 

FOR K a 1 TO IB: ZR = COS(IC * U): ZI = SIN(IC * U) 

FOR N a ic TO ID - 1: AR = Y1R(N + IC): AI a Y1I(N + IC) 

BR = ZR * Y1R(N + ID) - ZI * YlI(N + ID) 

BI = ZI * Y1R(N + ID) + ZR * Y1I(N + ID) 

Y2R(N) = AR + BR: Y2I(N) = AI + BI 

Y2R(N + NH) = AR - BR: Y2I(N + NH) = AI - BI: NEXT N 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2 : IB = 2 * IB 

FOR N = 0 TO NP - 1: Y1R(N) a Y2R(N) : Y1I(N) a Y2I(N): NEXT N: NEXT 

FOR N = NP TO 1 STEP -1: Y(N) = Y1R(N - 1) / NP: NEXT N: N1 = NH 

FOR L = 1 TO LP - 1: N1 a N1 / 2 

FOR N = 1 TO Nl: V(N) = .5 * (Y(2 * N - 1) - Y(2 * N) ) “ 2: NEXT N: 

VBAR a V(l) 

FOR N = 2 TO Nl: VBAR a VBAR + (V(N) - VBAR) / N: NEXT N 

FOR N = 1 TO Nl: Y(N) = .5 * (Y(2 * N - 1) + Y(2 * N) ) : NEXT N 

PRINT L, T* (2 * (L- 1)), VBAR: NEXT L: BEEP 


IB = 1 
L 

IB = 1 
L 

a 1 


Figure 31 . Program for random walk noise with x = 102.4 s and ALLFFT. 
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PROGRAM ALLFFT USING F 

N INT TIME 

.2 
.4 
.8 
1.6 

3.2 
6.4 
12.8 
25.6 

51.2 

10 102.4 


TO PERFORM CONVOLUTION 

VARIANCE 

394.7613 

610.4606 

1093.399 

1712.953 

3951.828 

12355.85 

15130.14 

38779.55 

64481.73 

117921.1 


Figure 32. Data for random walk noise with x = 102.4 s and ALLFFT. 
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CLS 

DIM X(2048) , X1R(2048) , X1I(2048), X2R(2048), X2I(2048), HM(2048), HP(2048) 
DIM Y(2048) , YlR( 2048) , Y1I(2048), Y2R(2048), Y2I(2048), V(2048) 

PRINT "PROGRAM ALLFFT1 USING FFT WITH SPECTRUM SHAPING": PRINT 

PRINT "N" , "INT TIME", "VARIANCE": PRINT : T = .2 

NP = 2048: LP - 11: NH « NP / 2: PI - 3.14159: U = 2 * PI / NP 

FOR N = 0 TO NH - 1: X(N) * 100 * (RND - .5): NEXT N 

FOR N = NH TO NP - 1: X(N) * 0: NEXT N 

FOR M » 1 TO NH - 1: HM(M) * 1 / M: HP(M) = 0: NEXT M 

HM( 0 ) « 0: HP(0) = 0: HM(NH) = 0: HP ( NH ) = 0 

FOR M = NH + 1 TO NP - 1: HM(M) = HM(NP - M) : HP(M) = 0: NEXT M 

FOR M - 0 TO NP - 1: X1R(M) = X(M): XII (M) = 0: NEXT M: IA » NP / 2: IB = 1 

FOR L = 1 TO LP: IC « 0: ID = IA 

FOR K = 1 TO IB: ZR = COS(IC * U) : ZI » -SIN(IC * U) 

FOR M > IC TO ID - 1 

AR = X1R(M + IC): AI = XII (M + IC) 

BR = ZR * X1R(M + ID) - ZI * X1I(M + ID) 

BI = ZI * X1R(M + ID) + ZR * XII (M + ID) 

X2R(M) = AR + BR: X2I(M) - AI + BI 

X2R(M + NH) = AR - BR: X2I(M + NH) = AI - BI : NEXT M 

IC = ID: ID * ID t IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR M ■ 0 TO NP - 1: X1R(M) = X2R(M): X1I(M) « X2I(M): NEXT M: NEXT L 

FOR M * 0 TO NP - 1: HR * HM(M) * COS(HP(M)): HI = HM(M) * SIN(HP(M)) 

Y1R(M) « X1R(M) * HR - XII (M) * HI 

Y1I(M) * X1R(M) * HI + XII (M) * HR: NEXT M: IA » NP / 2: IB = 1 

FOR L » 1 TO LP: IC * 0: ID = IA 

FOR K = 1 TO IB: ZR - COS(IC * U) : ZI ■ SIN(IC * U) 

FOR N - IC TO ID - 1: AR = Y1R(N + IC): AI = Y1I(N + IC) 

BR = ZR * Y1R(N + ID) - ZI * Y1I(N + ID) 

BI = ZI * Y1R(N + ID) + ZR * Y1I(N + ID) 

Y2R(N) = AR + BR: Y2I(N) = AI + BI 

Y2R(N + NH) = AR - BR: Y2I(N + NH) ■ AI - BI : NEXT N 

IC = ID: ID = ID + IA: NEXT K: IA = IA / 2: IB = 2 * IB 

FOR N - 0 TO NP - 1: Y1R(N) = Y2R(N) : Y1I(N) - Y2I(N): NEXT N: NEXT L 

FOR N - NP TO 1 STEP -1: Y(N) = Y1R(N - 1) / NP: NEXT N: N1 = NH 

FOR L ■ 1 TO LP - 1: N1 = N1 / 2 

FOR N - 1 TO Nl: V(N) = .5 * (Y(2 *N- 1) -Y(2 *N)) * 2: NEXT N: 

VBAR = V(l) 

FOR N = 2 TO Nl: VBAR = VBAR + (V(N) - VBAR) / N: NEXT N 

FOR N ■ 1 TO Nl: Y(N) » .5 * ( Y ( 2 * N - 1) + Y(2 * N)): NEXT N 

PRINT L, T * (2 ~ (L - 1)), VBAR: NEXT L: BEEP 


Figure 33. Program for random walk noise with x = 102.4 s and ALLFFT1 . 
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PROGRAM ALLFFT1 USING FFT WITH SPECTRUM SHAPING 
N INT TIME VARIANCE 


.2 

2 . 777584E-03 

.4 

5 . 458743E-03 

.8 

1 . 051909E-02 

1.6 

2 . 125172E-02 

3.2 

3.561426E-02 

6.4 

7 . 147162E-02 

12.8 

.1857618 

25.6 

.4147886 

51.2 

.3673056 

102.4 

.3560083 


Figure 34. Data for random walk noise with x = 1 02.4 s and ALLFFT 1 . 
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Figure 35. Plot for 
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